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ABSTRACT 

We introduce a novel implementation of orbit-based (or Schwarzschild) modeling that allows dark matter 
density profiles to be calculated non-parametrically in nearby galaxies. Our models require no assumptions to 
be made about velocity anisotropy or the dark matter profile. The technique can be applied to any dispersion- 
supported stellar system, and we demonstrate its use by studying the Local Group dwarf spheroidal (dSph) 
galaxy Draco. We use existing kinematic data at larger radii and also present 12 new radial velocities within 
the central 13 pc obtained with the VIRUS-W integral field spectrograph on the 2.7m telescope at McDonald 
Observatory. Our non-parametric Schwarzschild models find strong evidence that the dark matter profile in 
Draco is cuspy for 20 < r < 700 pc. The profile for r > 20 pc is well-fit by a power law with slope a = -1.0 ±0.2, 
consistent with predictions from Cold Dark Matter (CDM) simulations. Our models confirm that, despite its 
low baryon content relative to other dSphs, Draco lives in a massive halo. 

Subject headings: dark matter — galaxies: dwarf — galaxies: individual (Draco) — galaxies: kinematics and 
dynamics — Local Group 



1. introduction 

Understanding how dark matter is distributed in low-mass 
galaxies is central to the study of galaxy formation in the cold 
dark matter (CDM) paradigm. The first CDM simulations 
predicted that all dark matter halos share a universal density 
profile with a cuspy i nner slop e of a = dlnpoiw/dlnr = -1 
dNavarro et alj|1996bl hereafter iNFWf) . When observers be- 
gan studying low-mass galaxies, however, the y mostly found 
halos with a unif orm densi ty a = core (Burkert 1995; 
Persic et all Il996t iBorriello & Saluccl 1200 It Ide Blok et all 
200U iBlais-Ouellette et all 120011: ISimon et all 120051) . This 
disagreement between theorists and observers over the form 
of Pdm(j) became known as the core/cusp debate. 

Since the debate began, the number of profile parameter- 
izations used to describe low-mass galaxies by both theo- 
rists and observers has only increased. Observer s cham- 
pion empirical fits such as th e Burkert profile dBurkertl 
119951: ISalucci & Burkerj 120001) . cored isothermal models 
(Persi c et al.l 119961) or simply g eneric broken powe r laws 
koch et al.l boTfTHStrigari et alJ l2(M iWalker et al.l l2009h . 
Theorists have also introduced new fits to their simu- 
lated halos with varying, although still cuspy, inner slope s 
dNavarro et alJl2004t Istadel et alll2009t iNavarro et al.ll2010h . 
Modeling a galaxy with each of these parameterizations 
would not only be time consuming, but also biased if the true 
profile is unlike any of them. It is therefore preferable to use 
non-parametric methods to determine pDni(f)- 

Van den Bosch et al. (2006) first experimented with non- 
parametric orbit-based models by allowing the mass-to-light 
ratio M/L to vary with radius in the globular cluster M15. 
We introduce a similar modeling technique that uses axisym- 
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metric Schwarzschild modeling, combined with knowledge of 
the full line-of-sight velocity distribution (LOSVD) of stars, 
to break the well-known degeneracy between mass and orbital 
anisotropy. We demonstrate the capability of these models by 
applying them to the Local Group dwarf spheroidal (dSph) 
galaxy Draco. Draco is part of an interesting class of objects 
that are some of the most dark matter-dominated galaxies dis- 
covered. This makes differentiating between dark and lumi- 
nous mass in dSphs easier as the baryons have less of an ef- 
fect on the total density profile than they do in larger galax- 
ies. Recently , using improved data and modeling techniques, 
lAdams et al.l (120 121) found a cuspy dark matter profile in the 



low-mass galaxy NGC 2796 where previous studies found a 
core. Studies like these motivate us to investigate the dSphs 
with more sophisticated models. 

Our models represent a significant improvement over previ- 
ous work on dSphs as most studies use spherical Jeans mod- 
els (iGilmore et al1l2007l I Walker et al.ll2009llWolf et al.ll20l0h 
which require the modeler to make assumptions about the na- 
ture and degree of the anisotropy. These assumptions vary 
in complexity from simply assumin g isotropy, which can bias 
the estimate of a (IE vans et al.l 120091 ). to parameterizing the 
anisotropy as a g eneral function of radius (Strigar ret al.l2008t 
IWolf et al.|[20Tof) . Models that allow for more freedom in the 
anisotropy typically fall victim to the mass-anisotropy degen- 
eracy and cannot rob ustly determine the inner slope of pDM(r) 
(Walker et al. 2009). We hope to make a robust determina- 
tion of the inner slope in Draco with a suite of more general 
non-parametric Schwarzschild models. 

2. NON-PARAMETRIC SCHWARZSCHILD MODELS 

At the heart of our non-parametric technique is th e orbit- 
based modeling cod e developed by iGebhardt et al.1 (120001 
120031) . updated by iThomas et all (12004 120051). and de- 
scribed in detail in Siopis et al. (2009]). All orbit-based 
codes are base d on the princi p le of orbit superposition first 
introduced by [Schwarzschild dl979l). Similar axisymmet- 
ric co d es are used b y Rix et al.l d!997l)lvan der Marel et all 
d 19981). ICretton et all Id 19991) . and iValluri et al. (2004) while 
Ivan de n Bosch et al. (2008) present a fully triaxial modeling 
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code. The current Schwarzschild models that allow for dark 
matter do so by requiring the modeler to assume a parame- 
terization for the dark matter density profile Pdm(/)- Unfor- 
tunately, this parameterization is often exactly what we wish 
to determine. Current methods get around the circular logic 
of this dilemma by running models with different parameter- 
izations and comparing their relative goodness-of-fit with a 
X 2 test. Non-parametric modeling sidesteps the issue entirely, 
and lets the parameterization of Pdm(f) be output from the 
models, rather than input as a guess. 

The principle of orbit superposition works by choosing 
from a library of all possible stellar orbits only those orbits 
that best reproduce the observed kinematics of the galaxy be- 
ing modeled. If we know the mass density profile of the 
galaxy, and hence the potential, we can compute the appro- 
priate orbit library. However, since we do not know the po- 
tential of the galaxy, we must construct a number of mod- 
els with slightly different mass distributions and compare the 
goodness-of-fit of the resulting allowed orbits. The radial pro- 
file of the total (dark + stellar) mass density in a galaxy can 
be written as: 

M * ,s 

P(r)= — XL>(r)+p DM (r) (1) 

where M*/L is the mass-to-light ratio of the stars, v(r) is the 
stellar luminosity density profile, and Pdm(/) is the dark mat- 
ter density profile. In principle we know M^/L, which can 
vary as a function of radius, from stellar population models. 
We also know v{r) from the de-projection of the observed sur- 
face brightness profile. Our task is to construct orbit libraries 
for varying p{r) and match the allowed orbits to kinematics in 
the form of LOS VDs — the distribution of projected velocities 
observed. Some orbit libraries will contain orbits that do a 
good job at fitting the observed LOS VDs and others will not. 
The best-fitting model identifies the best-fitting p(r). Once 
we know this, we can invert Equation (Q~|i to solve for pDM(r). 
The trick is to vary p(r) in a systematic way. This is the 
principal difference between our new approach and standard 
Schwarzschild modeling which tries to vary p(r) by varying 
the parameters that define an assumed dark matter profile. 

To compute the orbit library for each model, we first cal- 
culate the potential. We assume axisymmetry and make 
use of the stellar ellipticity to define the density at angle 6 
in our meridional grid. The dark matter halo is assumed 
to have the same ellipticity as the stars. We solve Pois- 
son's equation for the potential associated with this density 
distribution by dec omposing p(r, 6) into spherical harmonics 
( Siopi s~et alj 120091) . With the potential known, we launch 
20,000-30,000 orbits and integrate their motion for roughly 
100 crossing times, storing position and velocity information 
at each timestep. 

Orbits in axisymmetric potentials respect three isolating in- 
tegrals of motion: energy E, the z-component of angular mo- 
mentum L z , and the non-classical third integral I3. By speci- 
fying all three of these quantities together, an orbit is uniquely 
defined. Unfortunately, there is no analytical form for /3 and 
it is generally not k nown a priori. We th erefore rely on the 
sampling scheme of Thomas et al. (2004) to construct an or- 
bit library which uniformly samples E, L z , and h and thereby 
contains all possible orbits for a given potential. 

Each orbit in the library is given a weight w,, and a set of 
w, are chosen so the observed kinematics are appropriately 
reproduced by the orbits which have been weighted, averaged, 



and projected. Quantitatively, we observe A^losvd LOSVDs 
in the galaxy at various positions. Each LOSVD contains A^ ve i 
velocity bins with uncertainties, so the number of observables 
the models must match to is given by the product A^losvd x 
iVyei- The goodness-of-fit of a model is judged by 

A'losvd Nvd / nobs _ pmod \ 2 

where tff and £fj oi are the / h velocity bin of the i th LOSVD 
from the observations and model respectively, and cr,y is the 
uncertainty in 

Given the freedom to choose from upwards of 10,000 or- 
bital weights to match only A^losvd x N ve \ ~ 100 observables, 
a standard \ 2 minimization routine can populate the distribu- 
tion function in any number of ways that introduce unwanted 
noise. To avoid distribution functions that are noisy or unre- 
alistic, however still consistent with the observables, we em- 
ploy a maximum entro p y smo othing technique d eveloped by 
iRichstone & Tremaind ([1988) and described in Siop is et all 
(2009J). Instead of minimizing x 2 , we maximize the objec- 
tive function 

^-^^[StJ-sX 2 (3) 

where N or b is the number of orbits in the library, and Af2, 
is the phase-space volume of the i th orbit. See Siop is et ail 
(2009) for a technical description of how we calculate phase- 
space volumes and maximize S. 

The first term in Equation <[3j is an entropy-like quantity 
and the second term is \ 2 from Equation (f2j). The parame- 
ter as controls which term influences S. For small a s , orbital 
weights are chosen to produce a smooth distribution function 
at the expense of reproducing the data. For large as, the data 
are well-fit by the model (x 2 is small), but the resulting dis- 
tribution function is likely not smooth. We determine the ap- 
propriate as for ea ch model using the scheme described in 
Siopis^eTarj d2009h . We start with as = and incrementally 
increase it until changes to \ 2 between successive iterations 
are small. Thus, the maximum entropy constraint serves to 
initialize the search for the minimum when a s = 0. By slowly 
increasing as, we drive down the importance of entropy to the 
fit until it no longer matters. 

2.1. Varying p(r) Between Models 

The major innovation of our new modeling technique is 
how we choose the density profile p(r) of each model. Current 
methods assume Pdm{f) and calculate p(r) from Equation ([1]), 
however this requires knowledge of the appropriate parame- 
terization for poM( r )- We use a fundamentally different strat- 
egy and divide p(r) into Abin discrete points whose value p at 
radius r, is labeled p,. The Abin points are spaced evenly in 
log r and connected by straight line segments. Each trial den- 
sity profile is now defined by the p, at each of the Ab; n bins. 
We run many models adjusting the values of the p,- so as to 
sample all possible density profiles. This strategy requires no 
assumptions to be made about the shape of p(r) or pDM(r), but 
it is computationally intensive for large Abin- 

The choice of A^in is a compromise between accuracy in 
reproducing p(r) and computational resources. Large values 
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FIG. 1 . — VIRUS-W IFU overlaid on top of an HST image from Segall et al. 1 2007). Red circles highlight fibers containing stars used in the determination of 
the central LOSVD. Note the HST PSF is significantly smaller than the typical 2" seeing at McDonald Observatory. 



of iVbin can make parameter space impossibly large, while 
small values can be overly restrictive on p(r). We have ex- 
perimented with iVbin= 5, 7, and 10. The added freedom with 
M>in=7 or 10 was not found to be worth the increase to the di- 
mensionality of parameter space. We have also tried connect- 
ing the pj with splines, but found the additional freedom pro- 
duced unrealistic density profiles. Concern over the smooth- 
ness of p(r) may be mitigated by the fact that p{r) only matters 
to our models in that it determines the potential. As the poten- 
tial is the integral of p(r), this introduces additional smooth- 
ness. 

We extrapolate the density at the outermost point as a power 
law with slope a^. The only parameters in the model are the 
Pi themselves and the extrapolation slope ctoo. The models 
also have the flexibility to add a central black hole of mass 
M, to the galaxy for future studies. 

2.2. Separating dark from stellar mass 

Once the best-fitting p(r) is found, the task remains still to 
recover the underlying dark matter density profile. This in- 
volves finding some other constraint on the stellar mass-to- 
light ratio. We can often determine M*/L from simple stellar 
population (SSP) models. The required input for SSP mod- 
els varies greatly, and different methods are appropriate de- 
pending on the galaxy modeled. For example, if spectra are 
available, stellar population synthesis models or Lick indices 
can be used. Lacking spectra, one ca n use the relations be - 
tween broad-band colors and M*/L dBell & de Jongl 120011) . 
In nearby dSph galaxies where individual stars are resolved, 
color-magnitude diagrams can be constructed to fit for age 
and metallicity with isochrones. We can also evaluate the ra- 
dial variation of M* /L as well without much additional effort. 
Spectral or photometric data need only be spatially binned 
with the same procedure repeated at each bin. Once M^/L is 



calculated, stellar density is simply the product of the (possi- 
bly radially varying) M*/L x v{r). 

3. APPLICATION TO DRACO 

We apply our new non-parametric Schwarzschild modeling 
technique to study the nearby Draco dSph. Draco is a satel- 
lite galaxy of the Mi lky Way orbiting at a d istance from the 
sun of only 71 kpc (lOdenkirch en et all2001l) . At this distance 
individual stars are resolved even with ground-based observa- 
tories. Consequently, the data we use are radial velocities of 
individual stars. Radial velocities are available for 158 stars in 
Draco, and we present radial velocities from new observations 
of an additional 12 stars near the center of Draco. 

We choose Draco because it is the most dark matter- 
dominated of the "classical" (pre-SDSS) dSphs. We can 
therefore differentiate between dark and luminous mass more 
easily since the baryons contribute less to the total density 
profile than they do in larger galaxies. Consequently, we can 
absorb larger uncertainties in M* /L. The primary science goal 
of this work, and a future study of all dSphs, is to determine 
the functional form of the dark matter profile in dSphs and 
compare results to theoretical predictions by CDM. 

3.1. Data 
3.1.1. Kinematics 

We use a combination of published radial velocities and 
new observations for kinematics in D raco. Data exist at larger 
radii for 158 stars (Kleyn a et al.120021) . but we wish to explore 
the central region of Draco in order to have the best constraint 
on the inner slope of pDiw(r)- 

To accomplish this, we observe the center of Draco 
with the VIRUS-W integral field unit (IFU) spectrograph 
(Fabri cius et al.ll2008l) on the 2.7m Harlan J. Smith telescope 
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FIG. 2. — Color-magnitude diagram of stars near the center of Draco. Col- 
ored asterisks are stars we observe, coded according to their offset from 
Draco's systemic velocity V n . s . Red stars have \V— < 30 kms~', blue 
stars have \V — V S y S \ > 50 kms~', and the green star has a radial velocity 
between 30 and 50 km s -1 of V„ s . 



at McDonald Observatory . This instrument allows for a 
high density of stars to be observed simultaneously, but with 
the drawback that fibers are not positionable. There are 267 
fibers that cover the 105" x 55" field of view with a 1/3 fill 
factor. We observed the spectral region covering 4900A to 
5500A with a resolving power R ~ 9000. 

The observations took place over the first half of 5 nights 
from 2011 August 1-5 in excellent conditions. Seeing was 
typically 2"or better, which is smaller than the 3 "2 diam- 
eter fibers. The standard battery of bias, Hg-Ne arc lamp, 
and twilight calibration frames were taken at the start of each 
night. We use an early implementation of the Cure data re- 
duction software. Cure is being developed as the pipeline 
for the Hobby-E berly Dark Energy Experiment (HETDEX) 
dHill et al.ll2006l) . We briefly describe steps taken to reduce 
the VIRUS-W data. A detailed description of Cure is beyond 
the scope of this paper. 

We perform standard CC D processing steps, usi ng the fit- 
stools package (described in lGossl & Riffeser 2002), to create 
master bias, twilight fiat, and arc lamp images for each night. 
We use twilight fiats in combination with arc lamp images 
to determine the distortion solution-a two-dimensional map 
which translates the (x,y) position of a pixel on the CCD to a 
fiber number and wavelength. 

Our science frames consist of 15-minute integrations of a 
single pointing of the central part of the galaxy. Prior to ob- 
serving, we determined the optimal position of the IFU by 
examining Hub ble Space Telescop e (HST) photometry of the 
central region ( Seg all et al.l 120071) . With accurate fiber and 
star positions, we determined a pointing that maximizes the 
number of bright stars on fibers (see Figure [TJ. There are 57 
science frames with this pointing, totaling roughly 14 hours 
of integration. 

We apply each night's distortion solution to the science 
frames yielding rectified, wavelength-calibrated frames. We 
then collapse and median-combine these science frames along 
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FIG. 3. — LOSVD generated from the discrete velocities of 29 stars 

with the twilight fiat frames. Each night's stacked science 
frame is divided by the appropriate master fiat for that night. 

Since the majority of the 267 fibers in the IFU are on empty 
sky, we are able to calculate an accurate sky model directly 
from each science frame. We compute this sky model for each 
fiber on the chip using a moving-window average of 20 nearby 
fibers. We subtract the sky model from each frame, and the 
resulting sky-subtracted frames for each night are median- 
combined. 

We extract 1-D spectra from 17 fibers containing stars. Star 
2 in our sample is used as a velocity standard since it is 
the brighte st member st ar with known radial velocity from 
lArmandroffeTal] 0995). We cross-correlate the other 16 
spectra to star 2 using the IRAF task FXCORR. By cross- 
correlating to the spectrum of a star with known heliocen- 
tric radial velocity in Draco, we automatically remove the 
contribution from Earth-Sun motion. We perform the cross- 
correlation analysis on the combined image, and in doing so 
introduce a small bias due to the change in the heliocentric ve- 
locity correction over the course of the observing run. How- 
ever, the magnitude of this change is only 0.1 km s~ , much 
smaller than our uncertainties. 

We list the heliocentric radial velocities and Tonry-Davis 
Rtd values determined for the 12 stars we report as mem- 
bers in Table Q] The Tonry-Davis value indicates the relative 
strength of the primary peak in the cross-correlation function 
to the average ( Tonr v~& Davisll9"79l) . The right ascension and 
declination given for each star in TableQ]indicate the position 
of the center of the VIRUS-W fiber containing that star. 

To determine membership f or th e 17 stars, we use 
the photometry of Sega ll et ail ( 120071) to produce a color- 
magnitude diagram (CMD). Figure [2] presents the result- 
ing CMD, where the colored symbols indicate observed 
stars. We also group the stars according to their offset from 
Draco's systemic velocity, which we assume is V sys = -293 
kms" 1 (Armandroff et al.1 ll995h . Stars with radial velocity 
offsets greater than 50 km s -1 are classified as non-members, 
while stars with offsets less than 30 km s -1 are categorized 
as members. The one star with radial velocity V - V sys = 
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TABLE 1 

Radial Velocities Obtained with VIRUS-W 



Star 


RA 


Dec 


Vhdio km s 1 


AV h(;1 i km s 1 




1 


17 h 20 m 14 s .76 


+57°54'32:'40 


-288.1 


2.57 


4.76 


2 


17 h 20 m 07 s .49 


+57°54'32"04 


-299. 1 1 


1.89 1 




3 


17 h 20 m 06\12 


+57°54'32"40 


-293.1 


3.99 


8.75 


4 


17 h 20 m 14 s .ll 


+57°54'23:'04 


-310.9 


3.35 


6.01 


5 


17 h 20 ra 12U0 


+57° 54' 13" 68 


-270.6 


3.37 


7.47 


6 


17 h 20 ra 16 s .78 


+57 54'59f'76 


-276.2 


1.91 


12.98 


7 


17 h 20 m 08U4 


+57°55W12 


-258.4 


3.89 


7.87 


8 


17 h 20 m 16 s .44 


+57 o 54'55f'08 


-293.2 


6.05 


8.01 


9 


17 h 20 m 07 s .80 


+57° 54' 55 "44 


-307.6 


4.51 


10.92 


to 


17 h 20 m 09 s .48 


+57°54'50f'76 


-277.7 


3.61 


8.02 


11 


17 h 20 m 19 s .10 


+57°54'46;'08 


-292.2 


3.23 


10.94 


12 


17 h 20 m 17 s .ll 


+57°54'46"08 


-277.8 


2.17 


14.39 



NOTE. — Heliocentric radial veloci 
VIRUS-W at the center of Draco 
1 From lArmandroff et alj fT995T) 

32.6 ±3.9 kms"'( green symbol in Figure [2]i is classified as 
a possible member. Possible and non-members are discarded 
from further analysis, leaving 12 member stars. Note that 
blind sigma-clipping retains these same 12 stars as members. 

We have individual radial velocities for stars at positions 
around the galaxy, but our models want the distribution of 
radial velocities at each position — the LOSVDs. We group 
the individual velocities into spatial bins and determine the 
LO SVD at each bin via an adaptive kern el density estima- 
tor (lSilvermanll 986; Gebhardt et al. 1996). In velocity space, 
this procedure replaces each of the N discrete observations 
with a kernel of width h and height proportional to N~ l h~ l . 
We use the Epanechnikov kernel (an inverted parabola) and 
sum the contribution from each discrete velocity to obtain 
a non-parametric representation of the LOSVD. The lcr un- 
certainties on the LOSVDs are calculated through bootstrap 
resamplings of the data (i.e. sampli ng with replacement 
from the N velocity measurements, see Gebh ardt et aljfl996l: 
iJOTdel&Gebhardt 2012). In Figure [3] we show an example 
LOSVD. 

We combine the new VIRUS-W data with 158 addi tional 
radial velocities from the literature dKlevna et al.ll2002h . We 
divide these 170 radial velocities into 8 radial bins of roughly 
20 stars each. LOSVDs are calculated for each of these 
bins, yielding kinematics coverage over the radial range 25"- 
1500" (8 pc-500pc). We fit Gauss-Hermite moments to the 
8 LOSVDs and plot the kinematics in Figure [4] This is only 
done for comparison purposes as the models fit directly to the 
LOSVDs. We compare the velocity dispersion as determined 
from the Gauss-Hermite fit with the standard devi ation of the 
indivi dual velocities (using the biweight scale; see Bee rs et all 
1990) in order to determine the best value for the smoothing 
width h. 

The issue of foreground contamination frequently comes 
up in the study of dSphs using individual radial velocities. 
There is always the possibility that some fraction of the ob- 
served stars are members of the Milky Way. These stars 
would be velocity outliers and therefore artificially increase 
the measured velocity dispersion or, in our case, the width of 
the LOSVD. Fortunately, the foreground Milky Way stars are 
well-separated in velocity space from the lKlevna et al.1 (120021) 
sample. Contaminants are also unlikely to have colors and 
magnitudes that place them on the red giant branch of Draco's 



for the 12 member stars observed with 



color-magnitude diagram. lLokas et al.1 ((2005) use these two 
constraints to estimate th at there are of order 1-2 Milky Way 
contaminants in the entire Klevna et al. (2002) data set. 

3.1.2. Photometry 

Our models are required to not only match the observed 
LOSVDs but also the three-dimensional luminosity density 
profile v{r). The first step in obtaining v{r) is to measure 
the two-dimens ional surf ace brightness profile. We use the 
photometry of Segalfetal] ( 120071) who publish a number den- 
sity profile of stars in Draco. This profile covers the radial 
range from 15"- 2400". We extrapolate the profile as a power 
law out to R = 6000" by fitting a constant slope to the profile 
in logarithmic space. To convert the number density profile 
to an effective surface brightness profile, we apply an arbi- 
trary zeropoint shift in log space until the luminosity obtained 
by integrating the surface brightness pr ofile is consistent with 
the observed luminosity (Mateol [19981) . We plot this surface 
brightness profile in Figure[5] 

We deproject the surface bright ness pr o file ac cording to 
the procedure detailed in Gebhardt et al. (1996). We as- 
sume surfaces of constant luminosity density v are coaxial 
spheroids and perform an A bel inversion. For Drac o we 
adopt an ellipticity of e = 0.3 (Odenkir chen et aLll2001l) . We 
assume an inclination of i = 90° for simplicity. Inclination 
is typically one of th e more difficult quantities to constrain 
(Thomas et al. 2007b). In addition to simplifying our models, 
assuming i = 90° provides the advantage that the deprojection 
is unique. For a detailed discussion of how uncertainties in 
vie wing angle and geom etry propagate through our models 
see lThomas et alj d2007al) . 

The resulting luminosity density profile we calculate has an 
average logarithmic slope (d In v jd In r) = -0.4 inside 50 pc. 
In Figure[5jwe plot v(r) and also illustrate the positions of our 
kinematics data. 

3.2. Models 

Our non-parametric models of Draco use M,i n =5 radial bins 
spaced equally in logr from 15"to 2000". We initialize our 
search for the minimum with a brute force method, construct- 
ing a coarse grid in dimensions from which we calcu- 
late all possible permutations of the AW parameters and the 
extrapolated slope a^. Additionally, we require the density 
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FIG. 4. — Gauss-Hermite moments fit to the 8 LOSVDs generated from 170 
radial velocities. The solid line is the result of our best-fit model. 

profile of each model to be monotonically decreasing or con- 
stant. This is a natural constraint, and it significantly lowers 
the number of models needed to sample parameter space. 

Once the models defining the coarse grid are evaluated, we 
employ an iterative sampling scheme to focus in on and de- 
fine the minimum in better detail. This method takes all the 
models with \ 2 within x\ m of the minimum \mm as starting 
points. For each starting point, a fractional step of size <5,- is 
taken above and below the initial value, one at a time, for 
all the density bins. If there is no change to Xmin' men ^ * s 
decreased. This procedure is repeated until <5; is less than a 
specified threshold. Additional models are also run as needed 
to fill in regions of parameter space that appear interesting. 

We do not attempt to fit for as we clearly do not have 
kinematics in that radial range to constrain the mass. Instead, 
we treat as a nuisance parameter and marginalize over it 
in our analysis. We restrict the value of the extrapolated slope 
to ckoo £ {-2,-3,-4} and every p(r) we sample has been run 
with each of these values. These slopes ar e representative of 
the isothermal, NFW, and Hernquist ( 1990) density profiles. 

Since Draco orbits within the dark matter halo of the Milky 
Way, it is probable that is has been tidally stripped at large 
radii. To account for this, the density is truncated at the tidal 
radius r, defined by 

For reasonable values of the Milky Way's mass M, Draco's 
mass m, and the Galactocentric radius of Draco's orbit D (as- 
sumed circular), Equation|4]gives an approximate tidal radius 
r, r; 3 kpc. We therefore truncate p(r) at this radius. We also 
assume the dark halo in Draco has the same flattening as the 
stars and therefore leave qoM fixed at 0.7. In the future we 
plan to investigate models with varying qoiw, however that is 
not the focus of this paper. 

4. RESULTS 

The x 2 curves for all the p, are plotted in Figure[6] Each dot 
represents a single model, and the red curve is a smoothed fit 




r (arcsec) 

FIG. 5. — Surface brightness profile E(r) (dashed) and deprojected luminos- 
ity density profile v{r) (solid) used in our models. Horizontal lines near the 
x-axis indicate the radial position of our kinematics bins. Numbers refer to 
the number of radial velocities used per bin. Note the central location of the 
new VIRUS-W data (innermost bin) in comparison to existing data. 

to the minimum. We obtain the red curve through a smooth- 
ing process that is similar to a boxcar average, except that we 
take the biweight of the 7 lowest x 2 values within the boxcar. 
This method is therefore less sensitive to outliers than a tra- 
ditional boxcar average. When determining a smoothed fit to 
the minimum, one is tempted to use only the points with the 
lowest x 2 - However, numerical noise causes models to scat- 
ter to both higher and lower x 2 m some bins. This is difficult 
to see by eye because scatter to higher values of x 2 causes 
the models to blend in with the black points in Figure [6] while 
scatter to lower x 2 makes models appear to stand out. The 
sliding biweight robustly picks out the center of this distribu- 
tion 

The red curve plots x 2 (pd f° r eac h radial bin, and therefore 
gives an indication of the model-preferred density at radius r, . 
We estimate the la uncertainties on each of the p, by deter- 
mining the portion of each parameter's x 2 curve, marginalized 
over all other parameters, that lies within A% 2 = 1 of the over- 
all minimum. Figure [6] shows this limit as a horizontal line 
whose intersection with the red curve indicates the la range 
of the density at bin i. In all further analysis we identify the 
midpoint of this range as the best-fitting value and report un- 
certainties as symmetric about this value. 

In two cases, bins 4 and 5, there are secondary minima that 
extend almost to the Ax 2 = 1 line but not quite. It is likely that 
with perfect coverage of parameter space the area between 
these minima would be filled in. However, available com- 
putational resources limit the extend to which we can sam- 
ple parameter space. In order to be more conservative in our 
analysis we fit a quadratic in log p to these minima, centered 
roughly on the midpoint between them (blue curves in Figure 
i>. 

The best-fitting model has unreduced X^ma = 9-1, and the 
number of observables our models fit to is v = A^losvd x N ve \ = 
8x15= 120. If we were to naively calculate a reduced x 2 , we 
would estimate x 2 = 0.08. This low value of x 2 results from 
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FIG. 6. — x 1 curves for all of the p\ parameters. Each black dot represents a single model (combination of py ,p2, . . . ps) and the red curve is a smoothed fit to 
the minimum. The red curve in any panel therefore is the x" curve marginalized over the other density points. The unit of density is Mq pc~ 3 . In panels 4 and 5, 
the blue curve is a parabola in log p that we use to interpolate between two local minima. 



an overestimation of the number of independent degrees of 
freedom v. The adaptive kernel density estimator we use to 
compute the LOSVDs introduces correlations among neigh- 
boring velocity bins, therefore reducing the number of truly 
independent degrees of freedom. 

To account for this, we consider the Gauss-Hermite pa- 
rameterizations of our best-fitting model (solid line in Fig- 
ure HI and input LOSVDs (points with error bars in Fig- 
ure HI). This model has x 2 u GH = 0-33 where vqu is 4 Gauss- 
Hermite parameters x 8 LOSVDs= 32. This x^ GH is still less 
than 1, however it is more consistent with previous studies 
dGebhardt et al.ll2003l) and may be du e to correlations amon g 
the Gauss-Hermite parameters (e.g. [Houghton et al. 2006). 
We use xl GH t0 ca l cu l ate tne appropriate scaling to apply to 
our models which use the LOSVDs in determining % 2 . We 
scale all un-reduced x 1 values by the factor xl /xl = 4-3 



4.1. Obtaining M^/L 

We have so far identified the best-fitting total density pro- 
file. In order to study the dark matter profile we must subtract 
the stellar density profile p*(r). This involves finding an in- 
dependent constraint on the stellar mass-to-light ratio M*/L. 
Using stars within the central 5' of Draco, we construct a 
g' — i' color-magnitud e diagram (CMD) f rom the photometry 
of lSegalletaT1d2007l) . We fit isochrones (iMarigo et al.ll2008l) 
to the CMD, corrected for Galactic extinction ( Schleg el et al.l 
[19981) . so that we may determine the age and metallicity of 
the stellar population. 

Figure [7] shows the CMD with our best isochrone fit. The 



red giant branch is well-defined, and we obtain a sensible fit 
with age f age = 12.7 Gyr a nd metallic i ty [Fe /Hl = -1.4. Us- 
ing the SSP models from iMarastonl ([2005) we are able to 
convert f age and [Fe /H] to a V-band stellar mass-to-light ra- 
tio M»/Lv = 2.9 ±0.6. Uncertainties in M*/Ly represent the 
spread in SSP predictions when different initial mass func- 
tions are assumed in the models. 

4.2. The Dark Matter Profile 

With M»/Ly determined from stellar population models, 
we can subtract p*(r) from the best-fitting total density profile 
obtained during the modeling procedure. We plot the result- 
ing dark matter profile in Figure [H] The red band is the 68% 
confidence band for each density point, marginalized over the 
others, and the gray band shows the 68% confidence band of 
all the parameters jointly (at A% 2 = 7.04). 

From Figure[8] it seems plausible that Pdm(/) can be fit by a 
power law of the form log pom = ol log r+f3 with the exception 
of perhaps the innermost data point. The slope of this fit a 
can be directly compared to both theoretical predictions and 
observations of similar dSphs. The innermost point, however, 
is puzzling. Its value indicates a large central density and a 
departure from the power-law nature of the outer profile. Fur- 
ther puzzling is that its point-wise uncertainty (plotted as a 
red error bar) indicates strong constraint despite the fact that 
we have no kinematic data in this region of the galaxy. We 
speculate that, in the absence of such data, models are able to 
arbitrarily increase the central density. Since the volume of 
this inner bin is small, the total amount of mass added is neg- 
ligible. With no kinematics in this region, models can easily 
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FIG. 7. — Color-magnitude diagram of stars within the central 5'of Draco. 
From left to right, we plot isochrones of (t age X 10 9 yr, [Fe/H]) = (11 .5,-1.6), 
(12.5,-1.4), and (13.5,-1.3). The solid red line is the (12.5,-1.4) isochrone 
we use when determining M»/Ly. 

absorb this mass without affecting x 2 ■ We therefore exclude 
the innermost point in all further analysis. 

The resulting power law fit to the outer four points is shown 
in blue in Figure [8] We characterize the uncertainty in this 
fit by constructing 1,000 Monte Carlo realizations with noise 
added to the density profile. We draw each point point i ran- 
domly from a Gaussian distribution with mean log p, and dis- 
persion given by the width of the lcr confidence band at point 
i in Figure [8] We repeat the fit for each realization, and de- 
termine the la uncertainties on a from the 68% span of this 
distribution. This procedure yields a = -1.0 ± 0.2. None of 
the 1,000 realizations has a slope a > -0.45 strongly indicat- 
ing that the galaxy is not cored for r > 20 pc. 

4.3. Orbit Structure 

Once we have determined the best-fitting model, we can 
calculate the internal (unprojected) moments of the distribu- 
tion function at each of the bins in our meridional grid. Of 
interest is the anisotropy in the velocity dispersion tensor, 
which we quantify with the ratio a r ja, — the ratio of radial to 
tangential anisotropy in the galaxy. We define the tangential 
anisotropy a t as 



a^^-iaj + al + vl) (5) 

in spherical polar coordinates where v^ is the rotational veloc- 
ity. Streaming motions in the r and directions are assumed 
to be zero. We plot (T r /<j, in Figure [9] Since the LOSVDs 
we use in Draco contain contributions from stars at all angles 
8, we average a r and all quantities in Equation (0 when cal- 
culating <r r /(T,. Consequently, we lose the ability to evaluate 
anisotropy as a function of 8. This can be avoided if bet- 
ter kinematics coverage is available, either through more stars 
with radial velocities in dSphs or two-dimensional integral- 



FlG. 8. — Best-fitting dark matter density profile in Draco. The red shaded 
region represents the point-wise 68% confidence band for Pdm(>") (Ax 2 = !)• 
with the solid black line derived from forcing symmetric logarithmic errors. 
The gray shaded region is the 68% confidence band on pDM(r) considering 
all parameters jointly (Ax 2 = 7.04). We plot the innermost point (excluded 
from all futher analysis) an an error bar with the same color scheme. The 
solid blue line is the best power law fit to the profile, and the dashed line 
shows an r _I NFW-like profile. We plot the best-fitting NFW halo from a 
small grid of parametric models as the dashed green line. Vertical lines along 
the x-axis indicate the radial range of our kinematic data. 

field spectroscopy in more distant galaxies. Fortunately most 
other large dSphs in the Local Group have many more radial 
velocities publicly available. 

We plot a r /a, in Figure [9] over the radial range that our 
LOSVDs sample. We determine the uncertainties in a,-/ a, 
by the maximum/minimum values of o>/ a t for models within 
A\ 2 = 7.04 of Xmin (If f°i" M>in+1 degrees of freedom). We 
find evidence for radial anisotropy at all radii, consistent with 
the "tidal sti rring" theory describing the origin of the Milky 
Way dSphs dLokas et alj|2010t iKazantzidis et al.ll201 ll) . Un- 
certainties are large on <r,-/a t , likely due to the small number 
of radial velocities available as kinematic constraint. To con- 
strain the anisotropy better, more radial velocities are needed. 

5. DISCUSSION 

5.1. Improvement over Parametric Methods 

Since we eventually fit our non-parametric dark matter pro- 
file with a power law, one can ask why we do not initially use 
a power law -parameterized profile. This would seem advanta- 
geous, especially given the large parameter space required by 
non-parametric methods. This reasoning, however, relies on 
the assumption that we know the profile is a power law a pri- 
ori. The point of this study is to relax this assumption and see 
what type of profile comes out of the modeling, rather than 
impose unjustified interpretation on the problem. It happens 
that Draco hosts a nearly power law density profile, but by not 
assuming this a priori we allow more general models to be ex- 
plored. As a rough check that our models have converged to 
a global minimum, we run a small grid of parametric models 
with an NFW dark matter density profile. The best-fitting of 
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FIG. 9. — Ratio of the radial to tangential components of the velocity dis- 
persion. Values of ay/ov different from unity indicate anisotropy. The black 
line is our best-fitting model. 

these models is plotted in green in Figure [8] 

5.2. Interpreting the Dark Matter Profile 

It is important to note that we only constrain the dark mat- 
ter density profile over little more than a decade in radius 
from 20-700 pc. One could easily imagine our power law 
fit changing from a = -1 to a core (a = 0) inside of r ~ 20 pc. 
Likewise, the slope may also change at larger radii than 
r ~ 700 pc without our knowledge. The NFW density pro- 
file has an outer slope a = -3 for r 3> r s , but our profile does 
not change slope within our model grid. This could indicate 
that r s 3> 700 pc, but without knowledge of the outer slope we 
cannot say with certainty that the profile is NFW-like. 

Recent cosmological A^-body simulations have been found 
to produce density profiles shallower than the tra ditional a = 
-1 cusps dStadel et alj|2009t iNavarro et al.ll2010h . Many au- 
thors suggest that dar k matter profiles are best parameterized 
by the Einasto p r ofile (INavarro et al.l2004tlMerritt et alJ2005t 
iGao et all 12008: Na varro et al.ll2010i) where the slope varies 
with radius according to a power law a(r) oc r" . These profiles 
can have shallower cusps than NFW, but do not have constant 
slopes over a large range in radius. Our non-parametric den- 
sity profile is well-fit by a single power law from 20 < r < 
700 pc, but, again, this is a fairly narrow range in radius. Our 
models cannot rule out an Einasto-like change in slope outside 
this radial range. More kinematics are needed to characterize 
the density profile at large and small radii. 

When calculating the potential, we allow the outer slope of 
p(r) to vary between 2 < Ooo < 4 for r > 700 pc, but, unsur- 
prisingly, we are unable to constrain a^. Tidal effects may 
also alter the shape of pDM(r) since Draco is orbiting within 
the dark matter halo of the Milky Way. The tidal radius cal- 
culated from Equation (0]i is sufficiently large that tides are 
unlikely to affect the stellar component, but pDM(r) at large 
radii could be affected. If this is the case, pDiw(r) would de- 
cline more steeply than expected and the total mass enclosed 



would be smaller than what we calculate. 

The cuspy a = -1 dark matter profile we find in 
Draco stands in contrast to many othe r observational stud- 
ies of dSphs that find a = cores (iGilmqre et al.1 120071 : 
iWalker & Penarrubial201 UUardel & Gebhardtl2012l) . The ef- 
fects of baryons are still not well-understood, and could po- 
tentially drive a to different values on a galaxy-by-galaxy 
basis. These effects are the sum of at least two competing 
processes. Adiabatic compression dBlumenthal et al.l fl986) 
draws in dark matter boosting the central pom and driving 
a to more negative values. On the other hand, feedback 
from star formation and supernovae can c ause strong outflows 
( INavarro et alj|1996al : iBinnev et al J 120011) which can in turn 
remove dark matter from the centers of galaxies, reshaping 
cuspy profiles into a = cores. 

In a recent paper, iGovernato et al.l (120121) use high resolu- 
tion cosmological A^-body simulations with a fully hydrody- 
namical treatment of baryons to test these two competing ef- 
fects in low-mass dwarf galaxies. They find that the cuspiness 
of the dark matter halo is directly related to the amount of star 
formation activity in the galaxy. This is expressed as a cor- 
relation between a and stellar mass M*. Their interpretation 
is that galaxies that form more stars (larger M„) have more 
supernovae and a greater potential to turn a cuspy dark matter 
profile into a core. Using their least-squares fit to the M*-a 
correlation, they predict a « -1 .3 (at 500 pc) for Draco's stel- 
lar mass. This is in approximate agreement with our measured 
value of a = -1. 

Perhaps owing to the lack of stellar velocities available in 
Draco compared to other dSphs, there are not many studies 
investigating its dark matter profile throu gh dynamical mod- 
els. A rough comparison can be made with Lokas et al. (2005) 
who fit profiles with an inner slope of a = -1 and an outer ex- 
ponential cutoff at large radii. They find a total mass-to-light 
ratio that varies with radius between 100 > M^/Ly > 1000 in 
the inner ~ 700 pc. These values are comparable to the total 
mass-to-lig ht ratio we ca l culate in the inner ~ 300 pc. How- 
ever, unlike lLokas etafl (120051) we do not impose an expo- 
nential cutoff in Pdm(f) at large radii. Our calculated M tot /Lv 
therefore rises sharply at large radii where the stellar luminos- 
ity profile is decreasing much faster than pDAi(r)- 

Importantly, M tol /L v > M*/L v = 2.9 ± 0.6 (the stellar 
mass-to-light ratio we derive from SSP models) at all radii. 
This means we can confidently state that Draco is dark matter- 
dominated at all radii, allowing us to easily absorb errors in 
M t /Lv from SSP models. In other words, when determin- 
ing pDiwir) from Equation (1) the uncertainty in p(r) dom- 
inates the uncertainty in stellar density since the product 
M^/L x v(r) is much smaller than p(r). This is one of the 
reasons we choose to test this non-parametric technique on 
Draco first. In the future we plan to extend this analysis to the 
remaining Local Group dSphs, which are also thought to be 
dark matter-dominated everywhere. 

5.3. Draco's Mass 

We plot the enclosed mass profile of our models in Figure 
[TOl The shaded region is the la confidence band derived from 
the extreme values of M(r) for all models within A\ 2 = 5.84 
of the minimum (lcr for Nu n = 5 free parameters, marginaliz- 
ing over aoo). The vertical ticks on the x-axis represent the 
radial extent of our kinematics coverage. From this plot it 
is apparent that, despite its low luminosity and stellar mass, 
Draco lives in a dark matter halo that is surprisingly massive. 
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FIG. 10. — (Top): Enclosed Mass profile of our best-fittin g model (black 
line) and la confidence region. The green point is the Wolf et al. (2010) 
mass estimator. (Bottom): Circular speed profile and lc confidence region. 
Colors are the same as above. Vertical tick marks on the x-axis represent the 
range of our kinematics coverage. 

An interesting comparison can be made with the brightest 
dSph Fornax, roughly two orders of magnitude higher in lu- 
minosity. If we compare the mass enclosed within a common 
physical radius of 300 pc, we find that for Draco M300 = M(r = 
300 pc) = 3.8!ftg x 10 7 M Q , and Uardel & Gebhardtl (120121) 
measure M300 = 3.5^'jj x 10 6 M Q for Fornax. Of course, For- 
nax is much more extended than Draco so it is sensible to 
also compare the mass enclosed within the deprojected half- 
light radius of each galaxy's stellar component. For Draco we 
measure M1/2 = M(r = r e ) = 1.6^| x 10 7 M Q , and in Fornax 
Uardel & Gebhardtl (12012I) measure M I/2 = 5.8^ x 1O 7 M . 
We would prefer to compare the total mass of each galaxy, 
but there are no kinematic tracers far enough out in the halo 
that the density profile declines sharply enough to keep mass 
finite for any dSph. Consequently, we cannot constrain the 
total mass observationally and we must rely on comparions to 
simulations (Section 5.4). 

We also use our dynamical models to compare our measure- 
ment of M\ n with the convenient mass estim ator proposed by 
Wolf et al. (2010} (see lWalker et al.l2009l and CappellarietajJ 
20061 for similar formulae). This formula relates M]/ 2 to the 
directly observable luminosity-weighted line-of-sight veloc- 
ity dispersion < crf m > and projected half-light radius R e . 
Tlie lWolfetaI1(l2010l) mass estimator is written as: 

M l/2 ^4G- l R e <a 2 LOS > (6) 

and lWolfetaHdlQloh give a theoretical argument for why the 
effect of anisotropy is minimized near r e for a variety of stellar 
systems in spherical symmetry. 

For a more fair comparison of Equation (|6]l to our mod- 
els we calculate M l / 2 from our data set, not the value listed 

in IWolf et all d20Tol) . We use < a\ os >= 1 1 .3 ± 1 .6 km s" 1 , 
calculated directly from our data in Figure [4] as well as 



R e = 158.1 pc and r e = 205.2 pc which we derive from the 
photometry in Figure [5] This calculation yields an estimated 
M1/2 = (1.9 ±0.5) x 10 7 M Q , in excellent agreement with the 
mass calculated from our models. We plot the estimated My/2 
as the green point in Figure [TO] 

5.4. Comparing Draco to CDM Simulations 

We can also gain insight into the properties of Draco's 
dark matter halo by examining the circular speed profile V c (r) 
plotted in the lower panel of Figure [10] . The green point 
plotted is Vi /9 = ./GMi n/c~ft = 20.0 ± 2.6 km s _1 using our 
value of the [Wolf et al.1 (120101) ma ss estimator. In a recent 
paper, Boylan-Kolch in et al.l (120 12h match the observed V\/2 
of Local Group dSphs to subhalos around a Milky W ay- 
like halo in the Aquarius Simulation (Springel et al. 2008) to 
derive constraints on each dSph's maximum circular speed 
Vmax — a quantity dir e ctly r elated to the total halo mass. 
iBoylan -Kolchi n et al.l (120121) find that this estimate of V max 
is usually 20 - 30 km s -1 smaller than the V max they obtain 
through abundance matching. These results lead them to con- 
clude that the Local Group dSphs are dynamically inconsis- 
tent with the types of halos they are predicted to inhabit from 
abundance matching. 

We are in a position to investigate this claim directly in 
Draco. We do not need to match our V\/2 to simulations in 
order to gain knowledge of V c (r); we calculate the latter di- 
rectly, and not just at the half-light radius. Interestingly, much 
of our circul ar speed profile lies above t he V max = 20.5^ g 
predicted by IBoylan-Kolchin et alj (120121) . At r = 500 pc, 
the radius where we run out of kinematic tracers and can 
therefore no longer robustly constrain the mass, we find V c = 
34.6!^ km s" 1 . We can take the lower bound of V c here as 
lower iimit on V max > 26.4. The scaling relations bet ween to- 
tal mass and V max for subhalos (Springel et al. 2008) imply a 
lower limit on Draco's total mass of M > 1 .0 x 10 9 M Q . 

Ours is not the first study to suggest that Draco l ives in a 
halo with such a large mass. Penarrubiaet al. (2008) demon- 
strate that a family of NFW halos with varying V max and r max 
are consistent with the stellar kinematics of any King model 
embedded in an NFW halo. They break this degeneracy by in- 
voking the correl ation between V max a nd r max found in CDM 
simulations (e.g. B ullock et al. 20011). Their study suggests 
that Draco is the most massive of the Milky Way dSphs with 
Vmax ~ 35 km s" 1 . 

The comparison between Draco and Fornax is interesting 
as the two galaxies are separated by almost two orders of 
magnitude in luminosity but may have similar masses. Since 
Draco's inner halo is nicely fit by the NFW density profile 
(Figure 0, we can rely on simulations to extrapolate a to- 
tal mass M > 1.0 x 10 9 M Q . However, multiple independent 
studies using differe nt methods sugges t that Fornax does not 
live in an NFW halo dGoerdt et all20 06: Wal ker & Pefiarrubid 
120111: Uardel & Gebhardtll2012l) . and we therefore should not 
use the NFW formalism to predict a total mass from its V max . 
Still, the similarity in the galaxies' values of Mj/ 2 and M300 
suggests that the simplest abundance matching models, which 
require a one-to-one mapping between luminosity and total 
mass, may not appropriately describe the dSphs. If Draco and 
Fornax do indeed have similar masses, despite vastly different 
baryonic properties, then there must be substantial stochastic - 
ity in the galaxy formation process at the dSph mass scale. 
Even without comparing to Fornax, it is clear that Draco's 
baryonic properties do not map in the expected way to its halo 
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